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Abstract 

Following the calculation of optimal energy transfer in thermal environment in our first paper 
(New J. Phys. 12, 105012), full quantum dynamics and leading-order 'classical' hopping kinetics are 
compared for the Fenna-Matthews-Olson (FMO) protein complex with the thermal bath described 
by classical white noise (the Haken-Strobl model) and by quantum Debye noise. Because the 
two-site quantum dynamics is fully captured by the hopping kinetics, the difference between the 
quantum dynamics and 'classical' kinetics is due to higher-order quantum corrections, which include 
non-local multi-site coherence, i.e. long-range coherence. We observe that higher-order quantum 
correction leads to small changes in the trapping time or in energy transfer efficiency (1-2% in the 
Haken-Strobl model and S 10% in the quantum Debye model), suggesting that quantum coherence 
does not extend beyond a pair of sites in FMO and the optimal efficiency in FMO can be explained 
by Forster rate theory. However, using the population flux network, we can identify significant 
differences in the major energy transfer pathway between hopping kinetics and quantum dynamics 
due to multi-site quantum coherence. (26% in the Haken-Strobl model and 32% in the quantum 
Debye model for the initial site at BChl 1). In addition, the golden-rule rate expression for the 
hopping rate provides a simple and reliable estimation on the stability of the energy transfer against 
the change of internal (Hamiltonian) and external (dephasing rate, trapping rate, etc.) parameters. 
The quantum-classical comparison with the removal of the bottleneck site, BChl 4, demonstrates 
quantum protection by the mechanism via multi-site quantum coherence. 



* E-mail: jianshu@mit.edu. 
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I. INTRODUCTION 



Natural photosynthesis is of particular interest due to its essential role as the energy 
source for life on earth. In the process of biological evolution over billions of years, photo- 
synthetic systems have developed optimal and robust strategies of converting solar energy 
to chemical energy. In the early stage of photosynthesis, solar energy is collected by pig- 
ments and transferred through light-harvesting protein complexes to the reaction center for 
the successive charge-separation. The energy conversion from photons to electrons is fast, 
robust, and nearly perfect in efficiency, although the total efficiency of photosynthesis is 
low. Understanding the mechanism of efficient energy transfer in natural light-harvesting 
systems can help developing low-cost and highly-efficient man-made solar energy apparatus, 
including photovoltaic devices and artificial photosynthesis [lj]. 

For a long time, energy transfer was considered as an incoherent process described by hop- 
ping kinetics with Forster rate constants. The Forster rate theory and its generalization have 
been a prevailing theoretical technique j^ ^]. In spite of the widespread success of the Forster 
rate approach, recent experimental advance has shown evidences of long-lived quantum co- 



aa 



herence in several natural light-harvesting systems, e.g., Fenna-Matthews-Olson (FMO) 
and phycocyanin 645 (PC645) [7]. A full quantum dynamic framework becomes necessary 
for studying coherent energy transfer. Many theoretical techniques have been developed to 
serve this purpose 29] . With temporal-spatial correlation for the protein environment, 



the generalized Bloch-Redfield (GBR) equation 8MU|, the hierarchic equation 



12 



and 



several other methods [20|, |27Jhave successfully predicted the long-lived quantum coherent 
phenomenon. Alternatively, the Haken-Strobl model and its generalization have attracted 
much attention due to its simplicity, although the bath noise is classical jol I231 26 , 28, 29 ]. 



Different theoretical methods have been tested in the simple two-site system [13|, |14| and 
other complex systems |30|, mainly focusing on the reliability of theoretical predictions. 

However, a systematic and comprehensive investigation is still needed to distinguish hop- 
ping kinetics and full quantum dynamics, with the goal of quantifying nontrivial quantum 
effects, i.e., long-range quantum coherence, in a general energy transfer network. In this pa- 
per, we will make the quantum-classical comparison in the FMO system with two different 
descriptions of baths: the classical white noise (Haken-Strobl model) and the quantum Debye 
noise. Our comparison is based on a self-consistent kinetic mapping of quantum dynamics, 
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and its detailed procedure was described in a review paper for the classical noise [23J and 



will be described for the quantum noise 



31] . Here we will use the leading-order kinetics: a 



hopping picture where the rate follows Fermi's golden rule. With a dipole-dipole interaction 
between two chromophores, Fermi's golden rule rate becomes Forster rate of energy transfer. 
In the standard fashion, such hopping kinetics is considered as a 'classical' description of 
energy transfer. In this paper, we will explore energy transfer in FMO using exact quantum 
dynamic equations and using Fermi's golden rule rate (i.e., Forster rate) to extract the dif- 
ference between full quantum dynamics and 'classical' hopping kinetics, which defines the 
contribution of nontrivial quantum effects: multiple-site coherence or long-range quantum 
coherence. 

r~| 

In the first paper of this series [9|, we applied the Haken-Strobl model and the GBR 
equation approach to optimize energy transfer with intermediate values for various bath 
parameters, such as reorganization energy, bath relaxation rate, temperature, and spatial 
correlation. In particular, we found an optimal temperature for efficient energy transfer 



in FMO. Our results have been verified by the hierarchic equation 32 



has been found in other conditions such as the spatial arrangement 



, and optimization 



chromophores 2 



and dipole orientation 



11] , site energy of 



To interpret the optimization behavior, 



we have proposed the concept of orthogonal (trapping-free) subspace, and determined the 



asymptotic scalings in the weak and strong dissipation limits j35j. Since two-site quantum 
coherence is included in Fermi's golden rule rate, 'classical' hopping kinetics can also predict 
the optimization behavior in many light-harvesting systems, as we will show in this paper. 
Therefore, our quantum-classical comparison in this paper is essential for identifying the 
contribution of nontrivial quantum effects to optimal energy transfer. Specifically, we will 
investigate two quantities: the trapping time and the population flux network. The former 
is directly related to the energy transfer efficiency as shown in our first paper 9|, whereas 
the latter is constructed by directional population flow for each two-site pair in the energy 



transfer network 
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A key advantage of natural photosynthesis compared to its artificial 
counterpart is the robustness against environmental variation and self-protection against 
damages. Here the quantum-classical comparison is combined with the stability analysis of 
energy transfer to quantify robustness in FMO. Our study thus provides a unique approach to 
understand the biological role of nontrivial quantum effects excluding the two-site coherence, 
different from other theoretical papers on time- dependent behaviors of quantum coherence 
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and entanglement 39l. 140 ] . 



The paper is organized as follows: In Sec. [TTJ we review the quantum dynamic framework 
for light-harvesting energy transfer. In Sees. IIIIHiVl we apply the Haken-Strobl model and 
compare the trapping times and the flux networks of FMO calculated from the leading-order 
hopping kinetics and full quantum dynamics. The sensitivity of parametric dependence for 
the trapping time is evaluated using 'classical' hopping kinetics. The robustness of energy 
transfer is explored by removing one donor site of FMO together with a quantum-classical 
comparison. In Sec. |V] we apply the Debye spectral density for the protein environment. The 
trapping time and the flux network are computed quantum mechanically using the hierarchic 
equation and classically using Fermi's golden rule rate (i.e., Forster rate). In Sec. I VI} we 
conclude and discuss our results of quantum-classical comparison and robustness analysis. 

II. LIOUVILLE DYNAMICS AND TRANSFER EFFICIENCY 

In this section, we review the theoretical framework of exciton dynamics, following the 
same notation as introduced in our first paper [9j. 

For each local chromophore (site) of the exciton system, a two-level truncation is reliable 
for the lowest electronic excitation. In consistence with low light absorption in natural light- 
harvesting systems, we consider the situation of single excitation, and then energy transfer 
dynamics can be studied in the subspace of single-excitation quantum states. Thus, we 



introduce a tight-binding Hamiltonian in the site ({|^)}) representation 41 j 



H = y% n |n)(n| + } j Jmn\m){n\, (1) 

where e m is the excitation energy at chromophore site m and J mn is the electronic coupling 
strength between the m-th and n-th sites. The system investigated in this paper is the Fenna- 
Matthews- Olson (FMO) protein complex with seven bacteriochlorophyll (BChl) sites ||], [f3, 



42M47||. and the Hamiltonian has been given in literature and can be found in our first 
paper [9]. 

For an exciton system, the time evolution of the reduced density matrix p(t) is governed 



by the Liouville equation [23|, |4l|, |48|, |49j 



pif) — ~ [£-sys + Arap + ^decay ^ ^dissp] pif) ■ (2) 
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The four Liouville superoperators C on the right-hand side of the above equation correspond 
to four distinct dynamic processes, which are discussed as follows. For an isolated system, the 
system Liouville superoperator £ sys is given by the commutator of the system Hamiltonian, 
AsysP = i [H, p] , and its explicit form in the Liouville space is 



sys\mn,kl 



(3) 



For conciseness, we neglect the reduced Planck constant h throughout this paper. 

The irreversible population depletion of the exciton system originates from exciton decay 
by the electron-hole recombination and energy trapping at the reaction center 0, The 
Liouville superoperators of these two processes are diagonal: [£deca y ]mn = kd- mn = {kd- m + 



k dtn )/2, and [C 



trapj ran 



k 



t:mn 



{h,m + k t ,n)/2, where kd, n and k t , n are phenomenological 



decay and trapping rate constants at site n, respectively. Here [C) mn = [C] m n,mn represents 
the diagonal element. In practice, we often assume a homogeneous decay process with 
kd-,n = kd- In the FMO system, BChl 3 is the trap site connecting to the reaction center, 
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28]. 



kt-n = k t 5 nj3 , and the trapping rate is set to be k t = 1 ps _1 

In addition to the above three dynamic processes, the excitation energy transfer is modu- 
lated by fluctuations due to the interaction between the exciton system and the protein envi- 
ronment. On the microscopic level, £dissp is evaluated using the explicit system-bath Hamil- 
tonian. Within this description, the linearly-coupled harmonic bath, Hsr = J2 n \ n )( n \B n , is 



widely applied, with B n the linear quantum operator of bath 
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491 ] . The dissipative dy- 



namics of system is then fully determined by the bath spectral density J(u). For simplicity, 
we ignore the spatial correlation of bath in this paper. Next we can apply quantum dynamic 
methods, e.g., the Redfield equation 50], the generalized Bloch-Redfield equation 10], 



and the Forster equation 0, |4|, under various approximations. For exponentially decay- 
ing bath time correlation functions, the hierarchic equation approach can provide a reliable 
prediction of quantum dynamics [12|-|l5|. 



Alternatively, we can view the s yste m-bath interaction as a time-dependent fluctuation 
on the system Hamiltonian 



21 



22 



511, i- e -, H(t) = H + 6H(t), with (5H(t)) = 0. The 



dissipative dynamics can be fully determined if all the time- averaged moments of SH(t) are 
resolved, which is usually an unfeasible task. In the extremely high temperature limit, 5H(t) 
behaves classically and the relevant second-order moment becomes real. One example of this 
approximation is the Haken-Strobl model where a classical white noise, (5e m (t)8e n (0)) = 
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r*5 m ,n£(t), is assumed on site energies |2l(. The dissipation Liouville superoperator becomes 
diagonal in the site representation, [£dissp] mn = (1 — <5m,n)r*, where T* is the pure dephasing 
rate. Since the Haken-Strobl model can be rigorously solved, it serves as the simplest model 
to examine our kinetic mapping of quantum dynamics. 

A key quantity of excitation energy transfer is the energy transfer efficiency q, which is 
the ratio of energy trapping at the reaction center. The mathematical definition of q is given 



by 
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28 



52| 



POO 

q= dt Tr {C tiap p(t)} = V] k t . n T n , (4) 
Jo 



where r n is the mean residence time at site n, r n = J °° dtp n (t). For an arbitrary vector 
X in Liouville space (e.g., the density matrix p), its trace is defined as Tr{X} = ^2 n X nn . 
In nature, spontaneous energy decay occurs on the time scale of nanosecond, much slower 
than the picosecond energy transfer process. The condition of kd ~ 1 ns _1 <C k t allows us 



to simplify the transfer efficiency to 9l. |23] 



where (t) = J2 n T n (k d = 0) is the mean first passage time to the trap state in the absence 
of decay (i.e., the average trapping time). The comparison of transfer efficiencies calculated 
from Eqs. (jl]) and 03) has been examined in our first paper [9J, and their excellent agreement 
over a broad range of T* proves the reliability of Eq. (jSJ). In this paper, we will ignore the 
energy decay process and focus on the average trapping time (t). Following the formal 
solution of Eq. (j2J), the average trapping time, 

(t)=Tr{£-V(0)} fcd=o , (6) 

is determined by the Liouville superoperator C = £ sys + £ tmp +£dissp and the initial condition 
p(0) = p(t = 0). For the FMO system, BChl 1 and BChl 6 connected to the baseplate are 
considered as two initial sites for energy transfer [15j. In our calculation, we consider two 
initial conditions at either BChl 1 (pi(0) = 1, initial condition I) or BChl 6 {pe(0) = 1, 
initial condition II). 
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III. QUANTUM-CLASSICAL COMPARISON IN THE HAKEN-STROBL 
MODEL 



A. Trapping Time 

In our first paper, we have demonstrated the genericity of optimal energy transfer by the 
competition of quantum coherence and bath-induced dissipation [9]. A remaining question 
is to identify contributions of nontrivial quantum effects. To do this, we systematically 
map the energy transfer process to a kinetic process. In this mapping, the leading order 
represents the 'classical' hopping behavior, and higher-order corrections represent nontrivial 
quantum effects. For the Haken-Strobl model, the kinetic mapping is solved in a recent 
feature article 23], following a stationary approximation for quantum coherence p mn (t). 



The same result is re-derived in the Laplace domain as shown in Appendix [A] The hopping 



rate between sites m and n is given by [23J, |53| 



k c - k c - 2Tmn I T I 2 (7) 

^mn ~ ^nm — p2 _i_ A 2 \° mn \ ' \' ) 

mn ' mn 

where the site energy difference is A mn = e m — e n and the overall dephasing rate is T mn = 
T* mn + k t - mn . In Appendix [S], we further prove that Eq. ([7]) can be recovered from Fermi's 



golden rule rate with a classical white noise. The resulting 'classical' rate equation reads 



23 1 



Pm — } (kmn-Prn ^nmPn) kt im P m , (8) 
nj^m 

where P m = p mm is the population at site m. Equation ^ can be organized into a matrix 
form as P = —{K c + K t )P, where [P] m — Pm is the population vector. The two rate 
matrices, K c and K t , are defined as \K c \ m ^ m) = -k^ n , [K c ] niTl = J2 m (^ n ) k m,m and 
[Kt] m ,n = 8 m ,nkt;n- The average trapping time and time-dependent population in the leading- 
order hopping kinetics can be calculated in a similar manner as defined for the quantum 
model introduced in the previous section, giving 

(t)c = Y, T n= Et(^ C + K t y l P{U)] n . (9) 

n n 

To distinguish the trapping time calculated by full quantum dynamics and by 'classical' 
hopping kinetics, we denote the quantum results as {t®, and the classical results as 



{t% , (t)c}- In Ref. [23J, we have proved that (t)Q and (t)c are the same for the two-site 



system in the Haken-Strobl model. Thus, the difference between these two trapping times 
arises from multi-site quantum coherence. 

In Fig. [TJ we plot the average quantum trapping time (£)q calculated by Eq. (J5J) and the 
classical counterpart (t)c by Eq. (jH]) as functions of the pure dephasing rate T* with the two 
initial conditions for the FMO system. We observe that the values of (t) computed from the 
two different methods are close. For example, the relative difference between (£)q and (t)c 
is 2% under the optimal condition of T* pt = 175 cm -1 for the initial population at BChl 1, 
while the difference becomes less than 1% under r* pt = 195 cm -1 for the initial population 
at BChl 6. Overall, the relative trapping time difference is always less than 10% for T* <1 30 



cm 1 . 



Our result shows that full quantum dynamics and hopping kinetics lead to similar behav- 
iors in the trapping time and energy transfer efficiency In the Haken-Strobl model, energy 
transfer in FMO is controlled by the downhill pathway from BChl 6 to BChl 3, which does 



not need long-range exchange assisted by multi-site quantum coherence [461 ] . More impor- 
tantly, the trapping time is the sum of residence time at all the sites so that the cancelation 
from different sites can reduce the quantum effect. Unlike the classical white noise, a quan- 
tum colored noise at a finite temperature complicates the quantum-classical comparison, 
which will be discussed in Sec. |V] 

B. Flux Network 

To illustrate the contribution of multi-site quantum coherence, we construct the flux 



network defined by directional population flows 
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. In the 'classical' kinetic model, the 



integrated population flux from site n to site m is defined as, 

F c =k c c _ k c c / ln \ 

mn n mn'n n nm m' V w / 

The equivalent definition of population flux F® n in quantum dynamics is less straightforward. 
In Appendix (Bj we use the quantum kinetic rate matrix K® to derive a rigorous definition, 

Fmn = 2 Im [^mriTnrjJ j (H) 

where T® m = J °° p nm (t)dt is the decay time of off-diagonal density matrix element p mn (quan- 
tum coherence). For the two-site system in the Haken-Strobl model, we can prove that the 
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above two equations lead to the same result, consistent with the quantum-classical com- 
parison of the trapping time. In the over-damped limit, i.e., T* ^> \J\, quantum dynamics 
reduces to 'classical' hopping kinetics, and the two definitions become identical for multi-site 
networks. Away from this limit, we can compare the results obtained with Eqs. ffTUl) and 
(JTTJ), and use the difference to define the higher-order contribution, e.g., multi-site quantum 
coherence in the Haken-Strobl model. Equations (TTOT) and (JTTJ) are valid for any quantum 
dynamic networks, and they will be tested with the Debye spectral density in Sec. PVT . 

Figure [2] presents relevant population fluxes (> 0.05) calculated by full quantum dy- 
namics (F® n ) and hopping kinetics (FfJ, with two initial conditions in FMO. For each 
initial condition, its respective optimal pure dephasing rate T* t is used. Both quantum and 
classical flux networks show two identical dominating pathways: 1 — > 2 — > 3 (path A) and 



6 —7- (5, 7) 4 —7- 3 (path B). This result is consistent with 2D electronic spectroscopy 46]. 
With initial condition II, path A is nearly negligible whereas with initial condition I, path 
B contributes significantly with F® 4 = 0.40 {F^ 4 = 0.56). As shown in the next section, 
the two-pathway structure in FMO helps maintain its high efficiency even if one donor is 
removed. 

Compared to results of the average trapping time, the quantum-classical difference of 
fluxes is much larger. For the flux F31 directly from BChl 1 to BChl 3, the weak electronic 
coupling J13 leads to a small value of F^{ = 0.053 in the hopping picture with initial condi- 
tion I. On the other hand, multi-site quantum coherence allows the long-range population 



exchange [49J, and this quantum tunneling effect enhances the quantum flux more than 
twice, F^ = 0.123. For the entire network, we introduce the weighted relative difference 
between classical and quantum fluxes as 

2 V , \F Q - F c I 
XF = ~Y \F Q + F c\ ■ (12) 

The result is xf — 26% for initial condition I and xf — 7% for initial condition II. Both 
values are much larger than those for the trapping time (< 2%). The comparison of popula- 
tion fluxes quantifies the contribution of multi-site quantum coherence in the Haken-Strobl 
model. 

The quantum-classical comparison of population flux also clarifies structures of various 
energy transfer pathways. For the first path of 1 — > 2 — > 3, our Hamiltonian shows that BChl 
2 is the energy barrier along the transfer pathway. As we discussed, the quantum tunneling 



effect is relevant, assisting energy transfer. For the second pathway of 6 — > (5, 7) — > 4 — > 3, 
the energy downhill structure allows quick energy transfer even in the 'classical' hopping 
picture, and quantum effect is less relevant. The weighted relative flux difference xf is 
much larger in initial condition I than that in initial condition II, consistent with energy 
structures of these pathways. 

IV. ROBUSTNESS OF ENERGY TRANSFER IN THE HAKEN-STROBL 
MODEL 

A. Stability Against the Variation of System/Bath Parameters 

The leading-order kinetic network provides a simple estimate of the parametric depen- 
dence of the average trapping time and can explain the robustness of energy transfer via the 
insensitivity to changes in dephasing rate, trapping rate, and system Hamiltonian. Within 
the kinetic network, estimation based on the magnitudes of the hopping rate and trapping 
rate yields the average trapping time on the order of picosecond time scale, (t) ~ 10 ps. 
In comparison with the average decay rate on the nanosecond time scale, kd = 1 ns -1 , the 
average trapping rate is much faster and therefore the efficiency is close to one, q ~ 1. The 
small ratio of the trapping time and decay time, kd(t) ~ 0.01, suggests that significant drop 
in energy transfer efficiency will result from a change of two-order's magnitude in the trap- 
ping time. According to Eq. ([7]), changes of two-order's magnitude in T* and k t or changes 
of one-order's magnitude in the J and A disorders are needed to produce this drop, based 
on simple estimation. The energy transfer in FMO can thus resist a dramatic change in 
dephasing rate, trapping rate, and system Hamilton. 

B. Robustness Against the Loss of a BChl Chromophore 

Over the course of evolution, the FMO complex has achieved robustness such that the 
high energy transfer efficiency is retained even if one or two chromophores do not function 
properly. We investigate this feature by taking out one non-trapping site off the FMO 
network and then calculating the average trapping time under the optimal condition. As 
shown in Table I, removal of BChl 4 causes a noticeable increase in the average trapping 
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Site Removed 


1 


2 


4 


5 


6 


7 


complete FMO 


(i)c(initial condition I) (ps) 




10.57 


15.32 


9.60 


10.61 


9.38 


10.30 


(i)<2 (initial condition I) (ps) 




10.32 


12.41 


9.17 


9.84 


9.13 


10.10 


{t)c (initial condition II) (ps) 


7.47 


7.61 


16.12 


7.98 




7.78 


8.72 


(t)Q (initial condition II) (ps) 


7.51 


7.72 


14.31 


7.80 




7.63 


8.66 



TABLE I: The quantum ((t)g) and classical ({t)c) trapping times of FMO with the removal of 
one non-trapping site under the respective optimal dephasing rate for two initial conditions in the 
Haken-Strobl model. As a comparison, the results for the complete seven-site model are shown in 
the last column. 

time while removal of any other site do not have major effects and some can even enhance 
the transfer efficiency. The flux distribution in Fig. [2] indicates that BChl 4 is the bottleneck 
of the dominant pathway 6 — > (5, 7) — > 4 — > 3, where all fluxes in this pathway converge 
to site 4 before arriving at the trap site. However, even with BChl 4 removed, the increase 
in the average trapping time or equivalently the decrease in efficiency is moderate. This 
phenomenon is due to the two major pathways in the network structure in FMO. Without 
multiple pathways in energy transfer networks, a linear-chain system exhibits much stronger 
response since the single pathway would be blocked after the removal of a donor. In addition, 
the 12-20% decrease in the trapping time (see Table I) demonstrates that the quantum effect 
can help the FMO system further resist the damage on BChl 4. Next we plot the relation 
of (t) - T* in Fig. [1] using both quantum dynamics and hopping kinetics for the FMO 
system with BChl 4 removed. Our results show that the quantum trapping time is always 
smaller than the classical result when BChl 4 is removed, verifying the protection of quantum 
coherence. 



V. GENERAL QUANTUM-CLASSICAL COMPARISON 



Our method of kinetic mapping for the Haken-Strobl model can be extended to a general 
quantum network, and its detailed derivation will be left in a forthcoming pape r |3l| . The 
leading-order hopping rate, however, is the same as Fermi's golden rule rate |49(j, which 



becomes the Forster rate if the dipole-dipo 
non-interacting blip approximation (NIBA) 



31 
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e interaction J mn is applied [2|. Following the 



551 ] . Fermi's golden rule rate is explicitly 
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written as 

i>oo 

kg^ n = 2 1 J mn | 2 Re / dte iAnmt e~ [9mn( - t)+kumnt] , (13) 
Jo 

and 

f°° J(u>) 

9mn{t) — 2(1 — c mn ) / duj — — [coth(/3o;/2)(l - coswt) + isinut] , (14) 
Jo u 

where = 1/fceT, J(u) is the bath spectral density, and c mn is the bath spatial correlation 
between sites m and n. The negative correlation c mri (^ m ) = — 1 is used in the spin-boson 



model 
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551 ]. and the zero spatial correlation c mn = 5 mi „ is considered in this paper. 



Consistent with our first paper [9j, we use the Debye spectral density, 

J H = M-^> (15) 

with 0(u;) the step function, A the reorganization energy, and D the Debye frequency. 
Without trapping, the ratio of forward and backward hopping rate constants satisfies the 
detailed balance condition, k^n/k^m = exp[— f3(e n — e m )}. The hopping rates from Eq. ( TIB"]) 
are then used to compute the trapping time and population fluxes in 'classical' kinetics. 



411 ]. the bath time correlation function C(t) of the 



Applying the Matsubara expansion 
Debye spectral density follows 

C(t) = Y,U-+*f^\ (16) 

3=0 

where the zeroth decay rate is vq = D and all the other decay rates are the Matsubara 
frequencies v 3 >\ = 2nj/(3. The coefficients /J and /j can be determined accordingly 15]. 
For the exponentially decaying bath, quantum dynamics can be rigorously solved by the 
hierarchic equation approach in principle 12l-ll6|. Here we use the explicit form shown in 



Ref. 



15], with the trapping Liouville superoperator £ t ra P included for both the reduced 



density matrix and auxiliary fields. The high-temperature approximation, coth(/3u;/2) rj 
2 //3a;, is applied so that all the Matsubara frequencies v 3 ~>\ are ignored in Eq. (JIB"]) . To 
be consistent, the same approximation is used in computing hopping rates. Our quantum 
computation is truncated upto the 10th hierarchic order. Due to instability of Liouville 
superoperators, is evaluated by the time integral of p(t) for A > 15 cm -1 . The resulting 
trapping time and population fluxes correspond to full quantum dynamics. 
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A. Trapping Time 



Here we study the quantum-classical comparison with the change of A by fixing T = 300 
K and l/D — 50 fs. The dependence on other parameters can be explored similarly. Figure[3] 
presents the results of the trapping time computed using both quantum dynamics ((£)q) and 
hopping kinetics ((t)c) with two distinct initial conditions. We observe that (£)q and (t)c 
are still close (j$ 10% for A ^ 30 cm -1 ) although their relative difference is larger than 
that in the Haken-Strobl model. Full quantum dynamics does not necessarily lead to a 
faster energy transfer process. With initial condition I, (t)q is smaller than (t)c for A > 34 
cm -1 ; with initial condition II, (t)g is larger than {t)c in the complete range of A. The 
latter behavior is consistent with a recent calculation of energy transfer rate in the two-site 
system [13]. 



To understand our quantum-classical comparison, we need to clarify higher-order correc- 
tions to the leading-order hopping kinetics in a quantum network at finite temperatures. 
(1) As shown in the Haken-Strobl model, the first quantum effect arises from multi-site 
quantum coherence, which facilitates the barrier-crossing energy transfer but does not con- 
tribute to the downhill pathway in a significant fraction. (2) The second quantum effect 
comes from system-bath entanglement, which is the signature of a non-Markovian bath. 
Without instantaneous bath relaxation, the Born approximation breaks down and the sys- 
tem is retained in its previous state, slowing down energy transfer. A systematic approach 
of including bath relaxation in electron transfer has been shown in Ref. 55|. (3) Other 
than the above two effects, a more subtle quantum-classical difference results from the 
steady-state population distribution. Without trapping, hopping kinetics always imposes 
P n (t — > oo) oc exp(— (3e n ). Instead, the quantum steady-state distribution is evaluated 
by p(t — > oo) oc Tre{exp(— (3H tot )}, which deviates from the Boltzmann distribution result. 
The steady-state population of the lowest energy trap site can be decreased, implying a lower 
probability of energy being trapped (energy transfer efficiency). The interplay of these three 
effects suggests that quantum energy transfer can be much more complicated, compared to 
the leading-order hopping kinetics. The identification of various higher-order effects for the 
trapping time and transfer efficiency will be studied in the context of the generalized non- 
interacting blip approximation (GNIBA) 31[ and the high-order multi-chromophore Forster 
resonance energy transfer (MC-FRET) theory 56]. 
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B. Flux Network 



For the physiological condition of A = 35 cm -1 15], we compare the quantum and 
classical flux networks, as shown in Fig. |H Similar to the study of the Haken-Strobl model, 
the flux network analysis clarifies the contribution of multi-site quantum coherence. With 
the barrier-crossing pathway under initial condition I, the quantum and classical trapping 
times are nearly the same but the detailed flux networks can be quite different: 6 — > (5, 7) — > 
4 — > 3 is the major path with a ratio of F% 4 = 75% in 'classical' hopping kinetics, whereas 
3 ^— 1 —7- 2 3 dominates in quantum dynamics with + F® 2 = 71%. The population 
flow from BChl 1 to BChl 3 differs by three times, F^/F^ = 3.06, and the overall quantum- 
classical flux difference is xf — 32%. The switch of the major energy transfer path results 
from the quantum tunneling effect of multi-site quantum coherence. With the downhill 
pathway under initial condition II, the flux network structure is the same for both hopping 
kinetics and quantum dynamics, with a small difference of xf = 7%. 



C. Robustness Against the Removal of BChl 4 

To complete our quantum-classical comparison for a general quantum dynamic network, 
we study the stability of FMO after the removal of BChl 4. As shown in Fig. [31 the change of 
the trapping time (t) is small enough to sustain a highly efficient energy transfer. Consistent 
with the Haken-Strobl model, quantum dynamics always leads to a smaller trapping time 
than hopping kinetics. This behavior can be interpreted by replacing the removal of BChl 
4 with an infinite energy barrier, e(BChl 4) = oo, so that the quantum tunneling effect is 
significantly important. As shown by the flux network analysis, BChl 4 is no longer the 
bottleneck site under initial condition I. The trapping time (£)q is unaffected by the removal 
of BChl 4 over a broad range of reorganization energy. The stability analysis thus reflects 
the multi-site quantum coherence feature of the original energy transfer network. 



VI. CONCLUSIONS AND DISCUSSION 



In this paper, we continue our investigation of efficient energy transfer in light-harvesting 
systems and compare the prediction of the trapping time and the population flux calculated 
by full quantum dynamics and by 'classical' hopping kinetics. The classical white noise (the 
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Haken-Strobl model) and the Debye spectral density are applied to simulate the protein 
environment. The quantum dynamics under the quantum Debye noise is solved by the 
hierarchic equation, which compares well to the GBR equation used in our first paper 9| 
but improves the accuracy. Relative to the rigorous results of quantum dynamics, hopping 
kinetics is calculated by Fermi's golden rule rate. In principle, full quantum dynamics can be 
mapped to an equivalent kinetic network of population transfer by a systematic expansion. 
We have extended this procedure from the Haken-Strobl model 123] to a general quantum 
dynamic network, which will be shown in a forthcoming paper [31]. In this mapping, the 
leading-order hopping rate is equivalent to Fermi's golden rule (the same as Forster rate for 
the dipole-dipole interaction) and is taken as the 'classical' hopping limit. Therefore, our 
quantum-classical comparison is capable of systematically illustrating nontrivial quantum 
effects using higher-order corrections. In the Haken-Strobl model, higher-order quantum 
effects originate purely from multi-site quantum coherence 23j. For a general bath model 
such as the Debye spectral density, there exist additional contributions from bath relaxation 
(non-Markovianity) |55[ and the finite temperature effect 31]. 

Our investigation in trapping time demonstrates that hopping kinetics compares well 
with full quantum dynamics and the Forster rate theory can reliably predict optimal energy 
transfer in FMO. Two initial sites, BChl 1 and BChl 6, are used in our study of FMO, 
reflecting the BChl arrangement in natural FMO. In both the Haken-Strobl model and the 
Debye spectral density, the quantum and classical trapping times are close, differing a few 
percentages over a broad parameter range around the optimal and physiological conditions. 
To achieve a highly efficient energy transfer, a rough energy funnel structure is naturally 
selected by FMO. In the local site basis, the energy difference (A) between two neighboring 
levels is several times larger than their electronic coupling (J). The study in our first paper 
suggests that an intermediate bath-induced dissipation (r), is necessary for the optimal 
energy transfer. A crude estimation on the kinetic expansion parameter, J 2 /(A 2 + T 2 ) < 1, 
indicates that quantum-classical difference can be treated as a small correction to the leading- 
order 'classical' kinetics in the overall dynamic behavior, i.e., the trapping time and the 
transfer efficiency. 

However, nontrivial quantum effects can be fundamentally important in the detailed 
behavior of the energy transfer process. A better measurement of nontrivial quantum effects 
involves off-diagonal coherence (p mn (^m)) of the reduced density matrix. In this paper, we 
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propose the population flux, which is defined using the time integral of p mn in full quantum 
dynamics. Through the flux network analysis, we are able to extract the two major energy 
transfer pathways, 1 ->■ 2 -> 3 (path A) and 6 ->■ (5, 7) -> 4 ->■ 3 (path B), in FMO. The 
energy transfer through path A crosses a barrier at BChl 2; path B is a downhill structure, 
thus becoming the dominant pathway in hopping kinetics. With the initial population 
at BChl 1 and the physiological condition of the Debye spectral density for the bath, the 
quantum tunneling effect switch the dominant pathway from path B to path A. The trapping 
time and flux network with the initial population at BChl 6 are much less affected due to its 
downhill structure. The quantum-classical comparison of the flux network thus characterizes 
multi-site quantum coherence for various network structures. 

Using the leading-order 'classical' kinetics, we present a simple estimation on the stability 
of energy transfer against the change of internal and external parameters. The time scale 
separation of energy trap and decay processes, kd(t) ~ 0.01, is a key factor for FMO. 
Based on the estimation of hopping rate, a noticeable change in the transfer efficiency 
requires a dramatic change in the trapping time, which needs one or two orders of magnitude 
difference in various parameters. To simulate a permanent damage in FMO, we remove the 
bottleneck site, BChl 4, and explore the modified trapping time in both quantum dynamics 
and 'classical' kinetics. We observe that the multiple pathways help FMO sustain a less 
dramatic change in the trapping time, and multi-site quantum coherence further accelerates 
energy transfer. 

Our analysis is based on a physically-motivated self-consistent kinetic mapping of quan- 
tum dynamics. The quantum coherence has been discussed in the framework of the long-lived 



quantum beat and entanglement [33, |40j. Complementary to these studies, we provide a 



quantitative measurement rather than a qualitative description for nontrivial higher-order 
quantum effects. Using the two-pathway FMO as our model system, we reveal the contri- 
bution of multi-site quantum coherence and its dependence on the pathway structure. Our 



approach can be easily applied to other light-harvesting systems [7j, [121, 1521 ■ Specifically, 
multi-site quantum coherence can lead to various phenomena, e.g., quantum interference 
between various energy transfer pathways, quantum phase modulation of a closed transfer 
.cop, and .ong-range energy exchange by quantum tunnehng Q. To systematically study 
these nontrivial quantum behaviors as well as the bath relaxation effect [55(, we need to 
develop a more detailed separation procedure based on our kinetic mapping technique [31 ]. 
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Our study is limited in the initial condition with pure populations, and a more general case 
with initial coherence will be extended in the future. For the FMO system, different Hamil- 
tonian models have been applied in theoretical and experimental studies 
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581 ] . and will 



lead to different results, e.g., the site energy of the additional eighth BChl can be optimized 
close to the experimental value [33]. 
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Appendix A: Quantum Kinetic Rate Matrix of the Haken-Strobl Model 

In this appendix, we will provide an alternative method to derive kinetic mapping in the 
Haken-Strobl model, and prove that the leading-order hopping rate can be recovered from 
Fermi's golden rule rate with a classical white noise. 

For a classical kinetic network satisfying the master equation in Eq. (jSJ), the integrated 
residence time is calculated by the Laplace transform 

P( z = Q) = [K + Kt}- 1 P(t = 0), (Al) 

where P(z) = J °° dte~ zt P(t) and r n = [P{z = 0)] n . For a quantum dynamic network, we 
will generate its kinetic mapping by the constraint of the same integrated residence time 
P(z = 0). The Liouville equation, p{t) = —Cp(t), is rewritten as 

pp(t) = —£ S ys;PCPc(t) — Arap;PPp(0 

Pc(t) = -C sys - C ppp{t) - [C sys -c ■^-'dissp;C ^trap; c}pc(t) (A2) 

where the indices P and C represent diagonal population elements and off-diagonal coher- 
ence elements of the density matrix in the local site basis. The reduced density matrix p is 
separated into two block elements, pp and pc- Each Liouville superoperator (£ sys , £ t ra P , and 
£diss P ), is also separated into the block- matrix form. In the Haken-Strobl model, the dissipa- 
tion Liouville superoperator is described the pure dephasing constant P*, i.e., £diss P; c = T*. 
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Applying the Laplace transform, we obtain a closed form for the population vector, 

Pp(z) = {z — C S y S -pc [Z + £ syS ;C + T* + £ t rap;c] 1 £>sys;CP + £trap ; p} Pp(t = 0), (A3) 

where pp(z) = J °° dte~ zt pp{t) is the population vector in the Laplace z-domain. To derive 
the above equation, we presume zero initial quantum coherence, pc(t = 0) = 0. With 
two identities, P = pp and K t = Arap;t> the integrated residence time vector from the full 
quantum dynamics is given by 

P(Z = 0) = {-C syS]PC [C sys - C + T* + A^C-r 1 £sys;CP + K t }~^ P(t = 0). (A4) 

Comparing Eq. flAlj) and Eq. (]A4j) . we obtain the quantum kinetic rate matrix as 

K Q = —C sys -pc [£ sys; C + T* + £ t rap;c] 1 £ S ys;CP> (A5) 

which includes the leading-order 'classical' hopping rates and higher-order corrections from 
multi-site quantum coherence. The effective quantum rate is given by k^ntem) = ~ [K®)mn- 
In the leading order, we ignore the off-diagonal elements of C sys -c, i.e., C sys -c ~ ^sys;C — * 
iA mn . With the explicit form of £ sys in Eq. (j3J), this simplification allows us to obtain the 
'classical' hopping rate k^ n fc£ n ) in Eq. ([7]). Next we can expand C syS] c in the order 
of £gyg. c and systematically determine the complete kinetic mapping of the Haken-Strobl 
model. 

In a forthcoming paper [3]]], we will demonstrate the mapping procedure for the general 
quantum dynamic network. The leading-order hopping rate is given by Fermi's golden rule 
rate in Eq. (ITBl . For a classical white noise described by J{oS) = (j3T* /2tt)u, we ignore the 
constant imaginary part of g(t) and arrive at 

? (i)»2 [°° d J^ T *^ )uJ j- (l-cos ut)=T*t (A6) 
Vo ^ $u 

and 

fc^= 2|J mn | 2 Re/ dre-^"^)^ 2] J mn | 2 - ^ , (A7) 

J mn 1" '- i r?in 

which recovers the result in the Haken-Strobl model. 
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Appendix B: Classical Flux and Quantum Flux 



For a classical kinetic network, the integrated population flux F^ n is denned by the net 
population flow from site n to site m, 

F c = uc c _ k c a mi \ 

mn mn n nm m" V / 

The quantum population flux F® n can be similarly defined by replacing the classical res- 
idence time and the hopping rate k^ n with the quantum residence time and the 
effective quantum rate k ( ^ in given in Eq. flA5j) for the Haken-Strobl model. Next we will 
prove that this definition of quantum flux is equivalent to Eq. (TTTj) . 

Applying the Laplace transform to the second equation of Eqs. ( 1A2I) and setting the 
Laplace variable to be zero {z = 0), we obtain 

- C sys . CP pp{0) = [£ sys; c + + £ trap;C ] p c {0), (B2) 

where the condition of zero initial quantum coherence (pc(t = 0) = 0) is used. The above 
equation be further rearranged, giving 

K^P(0)=C sys , P cp c (0) (B3) 

with the help of K® defined in Eq. (1A5I) . The n-th vector element on the both sides of this 
equation is written as 

^ ] \^mn T n ~~ ^nm T m\ = * ^ ] {Jnm T mn ~ Jrnn T nm) i (B4) 
m(^n) m(^n) 

where the coherence decay time r mn = J °° dtp mn (t) is introduced. Since the indices, m and 
n, are arbitrary in the above summation, we resolve the quantum population flux as 

^ran = kmn T n ~ ^nm T m = 2 I m [JmnT^ m ~\ , (B5) 

with two identities, J nm = J^ n and T® m = [t^J*. Although Eq. ( ]B5j) is derived in the 
Haken-Strobl model, this method can be applied to an arbitrary quantum dynamic network, 
giving the same result. Equation ( 1B5I) is thus the general definition of quantum flux. 
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Figure Captions 

Fig. [TJ The average trapping time (t) vs. the pure dephasing rate T* in the Haken-Strobl 
model for a) the initial population at BChl 1, and b) the initial population at BChl 6. The 
solid curves are calculated from full quantum dynamics using Eq. ([H]), whereas the dashed 
curves are calculated from the leading-order 'classical' kinetics using Eq. In each figure, 
the lower pair of curves correspond to the seven-site FMO model, while the upper pair 
corresponds to the six-site FMO model after the removal of BChl 4. 

Fig. HJ The flux networks of FMO under the optimal pure dephasing rate in the Haken-Strobl 
model for a) initial population at BChl 1 (r* t = 175 cm -1 ), and b) initial population at 
BChl 6 (r* t = 195 cm -1 ). For each population flux, the upper number is obtained from full 
quantum dynamics, whereas the lower number is obtained from the leading-order hopping 
kinetics. 

Fig. [3j The average trapping time (t) vs. the reorganization energy A of the Debye spectral 
density for a) the initial population at BChl 1, and b) the initial population at BChl 6. 
The other bath parameters are given in text. The solid curves are results of full quantum 
dynamics, whereas the dashed curves results of the leading-order hopping kinetics. In each 
figure, the lower pair of curves (black) correspond to the seven-site FMO model, while the 
upper pair (red) corresponds to the six-site FMO model with the removal of BChl 4. (colored 
online) 

Fig. HJ The flux networks of FMO under the physiological condition (A = 35 cm -1 , 1/D = 50 
fs, and T = 300 K) of the Debye spectral density for a) initial population at BChl 1, and 
b) initial population at BChl 6. For each population flux, the upper number is obtained 
using full quantum dynamics whereas the lower number is obtained using the leading-order 
'classical' hopping kinetics. 
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initial condition I r (cm ) initial condition II 

FIG. 1: The average trapping time (t) vs. the pure dephasing rate T* in the Haken-Strobl model 
for a) the initial population at BChl 1, and b) the initial population at BChl 6. The solid curves 
are calculated from full quantum dynamics using Eq. © , whereas the dashed curves are calculated 
from the leading-order 'classical' kinetics using Eq. Q. In each figure, the lower pair of curves 
correspond to the seven-site FMO model, while the upper pair corresponds to the six-site FMO 
model after the removal of BChl 4. 
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FIG. 2: The flux networks of FMO under the optimal pure dephasing rate in the Haken-Strobl 
model for a) initial population at BChl 1 (T* t = 175 cm -1 ), and b) initial population at BChl 
6 (r* pt = 195 cm -1 ). For each population flux, the upper number is obtained from full quantum 
dynamics, whereas the lower number is obtained from the leading-order hopping kinetics. 
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initial condition I ?i (cm ) initial condition II 

FIG. 3: The average trapping time (t) vs. the reorganization energy A of the Debye spectral density 
for a) the initial population at BChl 1, and b) the initial population at BChl 6. The other bath 
parameters are given in text. The solid curves are results of full quantum dynamics, whereas the 
dashed curves results of the leading-order hopping kinetics. In each figure, the lower pair of curves 
(black) correspond to the seven-site FMO model, while the upper pair (red) corresponds to the 
six-site FMO model with the removal of BChl 4. (colored online) 
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FIG. 4: The flux networks of FMO under the physiological condition (A = 35 cm -1 , 1/D = 50 fs, 
and T = 300 K) of the Debye spectral density for a) initial population at BChl 1, and b) initial 
population at BChl 6. For each population flux, the upper number is obtained using full quantum 
dynamics whereas the lower number is obtained using the leading-order 'classical' hopping kinetics. 
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